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ABSTRACT 


A  general  purpose  nonlinear  programming  method  and 
computer  code  are  presented.  The  method  is  basically  heuristic 
but  ^extremely  simple  and  reliable.  Interactive  operation  of  the 
code  and  its  performance  on  several  test  problems  is  described. 
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BACKGROUND 


Quite  commonly,  problems  in  the  statistical  and  engineer¬ 
ing  sciences  require  the  optimization  of  some  function  for  their 
solution.  Such  problems  include  those  of  resource  allocation, 
portfolio  selection,  curve  and  surface  fitting  (e.g. ,  linear 
or  nonlinear  regression),  approximation  theory,  signal  processing 
algorithms  (optimizing  a  model  fit) ,  and  optimal  control  (via 
discretization).  In  addition,  systems  of  nonlinear  equations  can 
be  solved  by  minimizing  some  function  of  the  magnitudes  of  the 
residuals.  As  for  the  statistical  sciences  in  general  we  have 
the  following  recent  comment  of  C.  R.  Rao:  "All  statistical  pro¬ 
cedures  are,  in  the  ultimate  analysis,  solutions  to  suitably 
formulated  optimization  problems"  [1], 

Throughout  this  report  we  shall  be  concerned  with 
optimizations  of  the  particular  form: 

minimize  f(x),  x  e  ft,  (1) 

where  ft  is  a  solid  subset  of  some  finite  dimensional  Euclidean 
space,  describable  by  finitely  many  contraints.  More  precisely, 
we  assume  that  ft  can  be  expressed  as 

ft  =  |xeRn:  ai_<xi£bi,  l£i<n,  c^ <g^  (x)£d^. ,  l£j£m|,  (2) 

with  the  understanding  the  ft  so  defined  has  a  non-empty  interior 
(is  "solid").  No  qualitative  assumptions  on  the  objective 
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function  f  in  (1)  or  on  the  constraint  functions  g^  in  (2),  such 
as  convexity,  differentiability,  etc. ,  are  made. 

Problems  of  the  form  of  (1)  are  known  as  mathematical 
programs.  They  have  been  the  subject  of  an  immense  amount  of 
study  over  the  lust  forty  years,  beginning  with  the  well-known 
special  case  of  linear  programs,  wherein  the  constraint  and 
objective  functions  are  linear  functions  of  their  variables.  As 
a  result  a  considerable  body  of  both  theory  and  computational 
algorithms  has  been  developed,  usually  under  convexity  and/or 
differentiability  assumptions.  The  books  [2-8]  provide  a 
representative  description  of  these  developments;  in  addition, 
there  are  specialized  journals  such  as  Journal  of  Optimization 
Theory  and  Applications  (since  1967)  and  Mathematical  Programming 
(since  1970 ) . 

There  is  a  basic  dichotomy  in  programming  algorithms; 
they  may  be  designed  to  converge  to  local  or  global  minima. 

Recall  that  a  local  minimum  of  f  is  a  point  x^eft  such  that 
f(x)_>f(Xg)  for  every  x  in  ft  that  is  sufficiently  near  to  xQ, 
while  a  global  minimum  is  indeed  a  solution  of  (1) .  Local  minima 
can  be  characterized  or  at  least  partially  identified  in  various 
ways,  depending  on  the  nature  of  the  problem.  A  prototypical 
necessary  condition  for  x^  to  be  a  local  minimum  of  a  differentiable 
program  (subject  to  a  "constraint  qualification"  on  the  geometry 


of  ft)  is  due  to  Kuhn  and  Tucker  [9]  and  might  be  considered  to  be  the 
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fundamental  theorem  of  mathematical  programming.  In  its  geometric 
form  it  states  in  essence  that  the  negative  gradient  of  f  at 
must  belong  to  the  cone  generated  by  the  gradients  at  xQ  of  the 
active  constraints  there.  Further,  under  various  generalized 
convexity  hypotheses,  e.g. ,  that  f  should  be  pseudoconvex  and  the 
g^  quasiconvex,  any  local  minimum  is  actually  a  global  minimum; 
in  such  a  case,  then,  the  Kuhn-Tucker  condition  characterizes  the 
solutions  of  (1).  This  remark  applies  a  fortiori  to  the  case 
where  all  functions  are  actually  convex. 

Such  characterizations  are  sometimes  used  to  provide 
stopping  criteria  for  numerical  algorithms.  For  example,  should 
a  local  minimum  occur  at  an  interior  point  of  ft,  as  would  certainly 
be  the  case  if  there  were  no  constraints,  the  K-T  condition 
reduces  to  the  vanishing  of  the  gradient  Vf  of  f  there,  and  in 
practice  an  algorithm  might  be  terminated  at  a  point  where  |  |Vf j  j 
becomes  sufficiently  small.  Of  course,  such  a  point  will  not 
necessarily  be  a  local  minimum  nor  even  close  to  one  in  the  absence 
of  convexity  conditions.  Similarly  the  popular  "method  of  feasible 
directions"  [3,10]  attempts  to  terminate  at  a  Kuhn-Tucker  point. 

In  the  important  special  case  where  all  constraint 
functions  are  linear,  so  that  ft  is  a  simplex,  the  constraint 
qualification  is  automatically  satisfied.  In  this  case  several 
effective  algorithms  exist  provided  that  V f  is  available.  These 
include  the  method  of  feasible  directions,  the  gradient  projection 
algorithm  [11],  and  an  approximation  method  wherein  f  is  replaced 
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by  its  first  order  Taylor  expansion;  this  approximation  yields  a  linear  program 
and  the  algorithm  generates  a  sequence  of  such  programs  (12] . 

In  this  linearly  constrained  case  there  are  very 
efficient  finite  algorithms  available  provided  the  objective 
function  f  is  either  linear  or  else  convex  and  quadratic. 

Failing  this,  but  with  Vf  available,  the  various  methods  mentioned 
in  the  preceding  paragraph  generate  sequences  of  linear  programs, 
systems  of  linear  equations,  and  line  searches  (one-dimensional 
optimizations)  to  be  solved.  Yet  other  methods  generate  sequences 
of  quadratic  programs  to  be  solved  by  constructing  quadratic 
approximations  to  f  [13,  14].  When  Vf  is  not  available  but 
can  be  assumed  to  exist  and  be  smooth,  a  powerful  but  rather  in¬ 
volved  algorithm  has  recently  been  described  in  [15]. 

For  the  general  (nonlinear ly)  constrained  mathematical 
program  there  are  adaptations  of  the  above  algorithms  obtained 
by  linearizing  the  contraints.  There  is  also  a  variety  of 
indirect  methods  which  may  operate  outside  the  feasible  set  ft. 
Examples  are  cutting  plane  algorithms  and  the  use  of  penalty 
functions  [2,  4].  The  former  generate  a  sequence  of  linear 
programs  over  a  decreasing  sequence  of  simplices  which  enclose  ft; 
the  possiblity  of  doing  so  depends  on  a  standard  trick  of 
utilizing  a  linear  objective  function  at  the  cost  of  introducing 
an  additional  variable  and  constraint.  The  latter  are  functions 
that  are  large  off  ft  and  which  are  added  to  f ;  a  sequence  of  these 
new  objective  functions  is  then  minimized  without  constraints.  All  these 


methods  involve  assorted  numerical  difficulties  as  well  as  the 
possiblity  of  termination  outside  ft.  Should  the  objective 
function  be  undefined  or  meaningless  off  ft, this  outcome  is 
naturally  untenable. 

In  the  last  decade  there  has  been  upsurge  of  interest 
in  the  global  optimization  problem  per  se;  see  [16,  17,  18]. 

This  problem  in  full  generality  is  very  different  from  and  more 
difficult  than  the  problem  of  local  optimization  just  discussed. 
There  is  relatively  little  theory  in  support  of  proposed  methods 
unless  fairly  strong  restrictions  are  placed  on  the  problems 
to  which  they  are  applied. 

The  methods  may  roughly  be  classified  as  deterministic 
or  probabilistic,  although  many  methods  have  both  features.  Most 
of  the  effort  has  gone  into  unconstrained  optimization  or  else 
into  problems  which  only  involve  bounds  on  the  variables. 

Of  these  the  probabilistic  methods  appear  in  general 
more  attractive,  in  part  because  of  the  lack  of  hypotheses  that 
must  be  imposed  on  the  problem  structure  and  in  part  because 
they  seem  less  sensitive  to  increases  in  dimension.  Another 
difference  concerns  the  manner  of  validation  of  particular 
algorithms:  with  probabilistic  methods  it  may  be  possible  to 
give  a  theoretical  a  priori  measure  of  tne  probability  of  "success", 
while  for  a  deterministic  algorithm  such  a  probability  can  at 
best  be  estimated  only  through  extensive  testing. 
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Among  the  several  popular  probabilistic  methods  we 
mention  specifically  only  two:  the  multistart  method  and  the 
sample  function  method.  The  former  is  essentially  an  iterative 
cycle  of  random  searches  throughout  £2  followed  by  application  of 
a  local  minimization  algorithm.  The  latter  hypothesizes  that  the 
objective  function  f  is  a  sample  function  of  a  known  stochastic 
process  on  £2.  The  statistics  of  the  process  then  permit  new 
points  in  £2  to  be  generated  according  to  the  information  provided 
by  the  values  of  f  at  preceding  points.  Further  statistical 
testing  can  be  employed  to  assess  the  likelihood  that  f  could 
indeed  be  a  sample  path  of  such  a  process.  Usually  the  Wiener 
process  is  assumed.  At  this  date  the  theory  and  computational 
requirements  seem  to  limit  this  method  to  small  (perhaps  only  1!) 
dimensional  problems. 

To  conclude  this  brief  background  we  note  that  a  general 
summary  of  computer  codes  for  mathematical  programming  that  have 
been  tested,  documented,  and  are  available  to  the  public  occurs 
in  [19].  Somewhat  earlier  there  appeared  a  collection  of  FORTRAN 
listings  of  optimization  codes,  along  with  brief  descriptions  of 
the  algorithms  and  their  operations  [20],  The  potential  user 
is  left  to  make  his  own  choice  as  to  which  method  will  best  serve 
his  purpose  and,  unfortunately,  the  various  codes  are  not  of 
uniformly  high  quality. 
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II. 


METHODOLOGICAL  AND  PROGRAMMING  DESIDERATA 


The  authors'  experience  with  numerous  "real-world" 
optimization  problems  has  persuaded  them  of  the  need  for  a  simple 
reliable  optimization  routine  that  could  be  interactively  operated. 
Such  a  program  should  in  particular  be  applicable  to  an  essentially 
arbitrary  objective  function,  about  which  little  or  no  prior 
information  may  be  available,  save  a  procedure  for  its  evaluation. 
This  situation  occurs,  for  example^  when  the  objective  function 
represents  the  output  of  some  "black  box"  operation,  perhaps 
from  a  process  optimization  or  simulation  model.  In  any  event 
the  nature  of  such  a  function  cannot  be  well  understood,  gradients 
are  unavailable  except  via  finite  difference  approximations, 
and  convexity  (occasionally  even  continuity)  is  in  doubt.  It  is 
also  often  the  case  that  the  feasible  region  must  be  strictly 
respected,  in  that  nonfeasible  points  have  no  significance  for 
the  problem. 

We  feel  that  the  nonavailability  of  analytic  gradients 
precludes  the  use  of  the  various  local  search  methods  mentioned 
in  Section  I.  While  some  codes  have  been  developed  that  utilize 
difference  approximations,  this  is  in  general  well  known  to  be  a 
numerically  risky  and  often  unstable  procedure.  There  are  the 
twin  perils  of  too  large  a  differencing  interval  (leading  to 
truncation  errors)  and  too  small  an  interval  (leading  to  round¬ 
off  and  cancellation  errors) .  The  situation  is  not  hopeless 


[15,  21,  22]  but  it  is  subtle  and  one  we  think  worth  avoiding, 
if  possible. 

The  requirements  of  simplicity,  robustness,  and  generality 
lead  us  to  a  combination  of  global  random  search  and  local  pattern 
search.  This  last  term  describes  a  class  of  methods  wherein 
the  objective  function  is  systematically  evaluated  at  several  points 
in  a  neighborhood  of  a  current  point,  the  resulting  function 
values  are  compared,  and  a  new  improved  point  is  selected.  A 
finite  number  of  failures  to  locate  an  improved  point  results  in 
termination.  These  methods,  and  their  termination  criteria,  are 
for  the  most  part  heuristic,  unsupported  by  significant  mathematical 
statements  of  optimality.  Nevertheless,  they  can  usually  be 
expected  to  perform  well,  if  suboptimally ;  they  will  rather 
quickly  produce  a  "good"  function  value,  although  extreme  accuracy 
may  require  many  further  evaluations,  and  premature  stalls  are 
possible.  Exact  solutions  to  specialized  problems  (such  as  an 
extreme  point  solution  of  a  concave  minimization)  will  not  normally 
be  determined  precisely. 

These  last  observations  lead  to  an  important  caveat: 
if  a  program  with  special  structure  is  confronted  there  will  most 
likely  exist  a  more  effective  specialized  algorithm.  Of  course 
tracking  down  such  a  method,  especially  in  reliable  coded  form,  may  not 
be  easy  (the  sources  in  [19,  20]  should  be  helpful  for  this 
purpose) .  Another  practical  point  is  this:  computer  time,  even 
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at  several  hundred  dollars/hour  is  still  cheaper  than  programmer/ 
analyst  time  at  $10-20/hour.  Even  if,  say,  a  routine  utilizing 
gradient  (and  possible  Hessian)  data  is  available,  these  deriva¬ 
tives  must  still  be  calculated  and  encoded.  The  savings  might 
typically  amount  to  a  few  tenths  of  a  second  of  execution  time, 
probably  not  a  cost  effective  trade-off.  An  important  possible 
exception  to  this  advice  occurs  when  the  purpose  of  the  optimizing 
code  is  to  serve  as  a  subroutine  of  some  large  program,  which  is 
to  process  data  in  real  time,  and  to  apply  some  control  law.  In 
such  a  case  speed  and  efficiency  are  likely  to  be  of  real  concern, 
and  specialized  software  may  be  required. 

A  review  of  several  popular  methods  of  pattern  search, 
including  the  complex  method  which  underlies  the  approach  described 
below,  is  given  in  Chapter  VII  of  [4]  and  in  [6,7],  We  also  take 
note  of  an  earlier  two-stage  heuristic  method  for  unconstrained 
optimization  announced  in  [23]  .  It  consists  of  a  single  non-ran¬ 
dom  global  search  of  pre-determined  length  from  a  given  initial 
point  followed  by  a  local  pattern  search.  The  latter  is  of  the 
Hooke-Jeeves  type  [6,7],  with  a  modification  to  account  for 
curvature  in  the  gradient  path.  No  restart  or  other  interactive 
features  are  offered. 


III.  THE  PROPOSED  METHOD 

The  following  method  meets  the  criteria  we  have  just 
imposed:  it  is  simple,  versatile,  robust,  and  can  be  operated 

interactively.  It  requires  as  inputs  the  absolute  minimum  of 
information  that  is  needed  to  specify  a  mathematical  program  of 
the  type  (1) ,  (2) ,  plus  a  few  parameters  that  essentially  deter¬ 

mine  how  much  work  the  program  will  do.  We  also  emphasize  that 
the  proposed  method  requries  no  linear  algebraic  subroutines, 
such  as  matrix  inversions  or  linear  equation  solvers,  nor  does 
it  construct  derivative  approximations.  Our  code  can  thus  serve 
as  an  "off  the  shelf"  package  for  optimization,  model  fitting,  and 
nonlinear  equation  solving.  Several  examples  illustrating  the 
program's  use  for  such  applications  will  be  given  in  Section  V. 

In  essence  the  program  utilizes  the  fundamental  idea 
of  Box  [24],  wherein  a  random  finite  set  of  feasible  points  is 
generated  and  successively  deformed  (based  on  information  gleaned 
from  function  evaluations)  until  it  collapses  around  a  local 
minimum,  or  else  (on  occasion)  becomes  mired  in  an  area  where  the 
function  is  essentially  constant.  Several  modifications  of  this 
basic  iteration  have  been  made  to  prevent  various  possible  traps 
and  infinite  loops  that  are  possible  when  the  feasible  set  fails 
to  be  convex  or  even  connected.  Further  modifications  have  been 
made  to  provide  several  starting  and  continuation  options  for  the 
users,  as  well  as  printing  and  termination  options. 
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To  describe  the  program  in  greater  detail  it  is  convenient  to 
treat  separately  its  three  major  aspects: 

1)  the  basic  iteration; 

2)  termination; 

3)  initialization. 

We  will  then  make  some  comments  about  the  continuation  and  print¬ 
ing  options 

Referring  to  the  notation  of  (1) ,  and  assuming  that 
ftCRn,  we  are  given  a  finite  subset  K  of  ft,  termed  a  "complex". 
Typically  the  size  of  K  (denoted  KPTS  in  the  program)  is  strictly 
between  n  and  2n;  unless  otherwise  instructed  our  program  sets 
KPTS  =  least  integer  >_  !.5n.  Three  points  of  K  are  now  singled 
out:  XLO,  XNXTHI,  XHI,  so  that 

f (XLO)  =  min{f (X) : XeK } 

<_  .  .  .  f  (XNXTHI)  < 

f  (XHI )  =  max{f (X) :XeK}. 

The  centroid  XC  of {XeK:  X^XHI)  is  now  computed  and,  for  the 
moment,  is  assumed  to  be  feasible  (i.e.,  XC  eft).  Next  the  point 
XHI  is  reflected  through  XC  to  a  point  XNU  so  that  dist  (XNU,  XC) 
=  STEP* dist (XC,  XHI);  here  STEP  is  a  program  parameter  typically 
in  the  range  (1,2].  Thus 

XNU  =  XHI  +  ( 1+STEP )  (XC-XHI).  (3) 
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Now  if  it  happens  that  XNU  is  both  feasible  and  satisfies  f{XNU) 

<  f(XNXTHI)  we  define  a  new  complex  by  replacing  in  K  the  point 
XHI  by  XNU;  in  this  case  the  basic  iteration  is  complete. 

Should  either  of  the  two  requirements  on  XNU  fail,  XNU 
is  backed  halfway  towards  the  centroid  XC  along  the  line  through 
XC  and  XHI.  If  this  point  meets  both  requirements  we  substitute 
it  for  XHI  and  obtain  a  new  complex.  If  not,  we  halve  again. 

If  after  a  finite  number  (typically  5-10)  of  such  attempts  we 
are  still  unsucessful,  we  set  XNU  =  XC  and  commence  moving  XNU 
halfway  towards  XLO,  again  trying  to  meet  the  twin  requirements  of 
feasibility  and  lower  f-value.  If  a  finite  number  (typically  10- 
20)  of  such  attempts  are  still  unsuccessful,  the  algorithm  is 
terminated  if  the  last  move  resulted  in  feasibility  only.  If  not, 
we  reflect  the  latest  XNU  an  equal  distance  through  XLO  and  if 
both  requirements  still  cannot  be  met,  the  algorithm  is  again 
terminated.  In  all  these  cases  of  premature  termination  a  warning 
message  is  printed  indicating  either  an  infeasible  direction  (and 
suggesting  that  the  best  point  XLO  is  located  in  a  very  "thin" 
region  of  ft) ,  or  else  the  inability  to  decrease  the  function  value 
sufficiently,  in  which  case  the  algorithm  is  simply  "stuck". 

There  is  one  other  conceivable  difficulty  that  can 
arise  during  a  basic  iteration,  and  that  occurs  when  the  centroid 
XC  is  not  feasible.  Of  course,  this  could  only  happen  when  ft  is 
not  convex.  In  this  case  we  again  set  XNU  =  XC  and  commence 
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moving  XNU  towards  XLO  as  just  described,  repeating  exactly  the 
above  steps.  Hence,  either  a  new  complex  is  obtained  or  else 
termination  ensues  with  the  appropriate  warning  message. 

Earlier  variants  of  this  basic  iteration  are  described 
in  the  original  source  [24]  and  in  [4,  20].  A  flow  chart  of  our 
version  appears  in  Fig.  1  (wherein  the  identification  Point ■* — ►XNU 
should  be  made) . 

In  most  "reasonable"  optimization  problems  the  basic 
iteration  as  just  described  will  result  in  a  new  complex.  Con¬ 
sequently  we  can  expect  to  generate  a  sequence  of  complexes  and  so 
must  decide  when  to  cease.  Roughly  we  should  do  so  when  no  further 
progress  is  being  achieved.  In  practice  we  evaluate  this  con¬ 
dition  by  means  of  the  absolute  and  relative  oscillation  of  f  on 
the  complexes.  If  either  of  these  quantities  drops  under  given 
thresholds  for  a  finite  number  (typically  5)  of  complexes  in 
succession,  the  algorithm  terminates.  Otherwise  a  new  basic 
iteration  is  carried  out.  Of  course  an  upper  bound  on  the  number 
of  iterations  must  also  be  supplied. 

Thus,  assuming  that  new  complexes  can  be  generated  we 
control  termination  with  four  program  parameters:  ABSTOL,  RELTOL, 
NTOL,  and  NLOOPS.  If  for  a  sequence  of  length  NTOL  of  complexes 
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The  basic  iteration 


Ki  +  1'  Ki+2 ' * ‘ '  Ki+NTOL  We  flnd 


f(XHI)  -  f(XLO)  <  ABSTOL, 


or  failing  this,  if 


f(XHI)  -  f(XLO)  <_  RELTOL*  |  f  (XHI)  |  , 


(4) 


(5) 


the  program  terminates  with  a  convergence  assertion.  Otherwise, 
it  terminates  with  a  warning  when  the  number  of  basic  iterations 
exceeds  NLOOPS.  These  termination  procedures  are  displayed  in  Fig.  2. 

To  initialize  the  optimization  process  we  must  construct 
a  first  complex.  Our  program  permits  this  to  be  done  in  three 
ways : 

a)  by  direct  assignment; 

b)  randomly  about  a  given  feasible  point; 

c)  choosing  the  best  KPTS  points  from  a  preliminary 
random  search  of  ft. 

Option  a)  might  be  used  when  the  geometric  structure  of  ft  was 
known  and  simple  (e.g.,  ft=  {  (x^ ,  . .  .  ,  xn>  :  0<x^  ,  1)  ,  and  no 

information  about  the  minimum  of  f  over  ft  was  available.  Option 
b)  is  valuable  when  there  is  a  priori  information  about  f  to  the 
extent  that  we  believe  that  for  some  x^eft,  f(x^)  is  close  to 
min  f(ft)  .  The  remaining  KPTS-1  points  for  the  initial  complex 
are  then  produced  sequentially  by  first  drawing  a  point  randomly 
within  the  region  defined  by  the  explicit  constraints  {x:  a^£x^ 

l£i£n}  in  (2),  and  then  retracting  it  by  halves  towards  x ^ 
until  it  becomes  feasible.  If  a  finite  number  (typically  10-20) 
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Pig.  2.  Algorithm  termination 


of  these  retractions  fails  to  result  in  feasibility,  we  begin 
again  with  a  new  random  point. 

In  the  remaining  circumstances  we  are  effectively  ignorant  of  both 
the  geometry  of  fi  and  the  behavior  of  f.  If  so,  we  utilize  option 
c)  by  making  an  initial  random  search  of  fi  and  retaining  the  best 
KPTS  points  to  form  the  initial  complex.  The  only  complication 
here  is  the  number  of  random  points  (denoted  NRANPT  in  the  program) 
to  utilize.  Clearly  we  must  have  NRANPT  >_  KPTS,  bu".  also  in 
practice  we  must  confront  the  trade-off  between  the  cost  in  terms 
of  function  evaluations  of  taking  NRANPT  large  vs.  the  cost  in 
terms  of  low  probability  of  being  close  to  the  minimum  if  NRANPT 
is  small. 


One  way  to  sensibly  determine  a  value  of  NRANPT  is  to 
require  a  certain  probability  of  placing  at  least  k  points  inside 
a  small  neighborhood  V  of  the  global  minimum  of  f  in  Q.  If  the 
relative  volume  of  V  in  is  a,  then 

Prob  ( _>  k  points  in  V  in  N  tries) 


=  1 


=  1 


>  1 


N-i 


5  (*)‘* 

,,  ,  N  k-1  /  a  \ 

(1-a)  N  exp (lT^J 


N! 

(N-i) J 


so  that  to  make  thi-s  probability  exceed  3  it  suffices  to  choose 
N  so  large  that 
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...  .  N  .,k-l  a  _  ,  „ 

(1-a)  N  exp  1  1_3* 

In  the  special  case  where  k=l ,  the  "exp"  term  can  be  dropped 
from  the  left  hand  side. 

To  illustrate,  let  us  take  V  be  to  a  cube  of  side  length 
2r  centered  at  the  global  minimum.  In  general  V  need  not  exist 
if  the  minimum  is  too  near  the  boundary  of  ft.  However,  if  there 
are  no  implicit  constraints,  or  if  we  can  be  sure  for  other 
reasons  that  the  minimum  is  well  inside  ft,  then  it  makes  sense 
to  discuss  V.  Now  a= (2r) n/vol (ft )  and  if  k=l  we  need  to  choose 
N  so  that 


N  >  >  isaU-a.  .  T_ vol(n). 

■  log  (1-a)  —  —  ot  (2r)n 

The  least  such  N  are  displayed  in  Fig.  3  for  nominal  values  of  B 
and  n,  with  r  set  =0.1  and  vol(ft)=l.  One  can  easily  note  a  very 
prominent  "curse  of  dimensionality"  here.  This  concludes  our 
discussion  of  the  initialization  step  in  the  program. 

The  progress  of  the  algorithm  when  applied  to  a  par¬ 
ticular  problem  can  be  monitored  as  desired  through  the  use  of 
four  print  options.  These  are  set  by  assigning  to  a  variable 
IPRINT  values  =  0,1, 2, 3.  The  0  value  supresses  all  printing; 
this  might  be  appropriate  when  the  optimization  is  serving  as 
a  subroutine  to  a  larger  program.  The  1  value  causes  only  the 
final  answer  (best  point  and  least  function  value)  and  machine 
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1 


TR-576 ( 3 ) 


2 

5 

10 

ID 

• 

O 

18 

2166 

6,769,016 

0.9 

58 

7196 

22,486,183 

0.95 

75 

9362 

29,255,198 

Fig.  3.  Least  number  of  random  samples  to 
give  probability  B  of  at  least  one  sample 
falling  in  cube  of  radius  0.1  within  a 
unit  n-volume. 

execution  time  to  be  displayed.  The  2  value  results  in  display  of 
the  best  initial  complex  point  and  corresponding  function  value 
followed  by  each  new  improved  point,  the  corresponding  lower  func¬ 
tion  value,  and  the  loop  number  w  en  this  improvement  occurred. 
Finally,  the  3  value  causes  the  entire  initial  complex  to  be  dis¬ 
played,  followed  by  further  information  for  every  iteration, 
whether  or  not  an  improvement  was  achieved.  Specifically,  the  new 
point  and  its  function  value  are  displayed  along  with  the  index  of 
the  point  in  the  previous  complex  which  it  replaces,  and  the 
number  of  function  evaluations  required  to  obtain  the  new  point. 

It  remains  to  discuss  the  continuation  options.  Recall 
that  the  program  will  terminate  for  one  of  three  reasons:  con¬ 
vergence  criterion  met,  number  of  iterations  loops  exceeded,  or 
prematurely,  because  of  failure  to  construct  a  new  complex.  Use 
of  the  1PRINT  =2  or  3  option  will  have  enabled  the  user  to  gain 
an  impression  of  the  progress  made  and  its  rate  since  initializa¬ 
tion.  If  this  progress  is  satisfactory  but  the  convergence  criterion 
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is  unmet,  the  user  will  probably  wish  to  continue  from  the  most 
recent  complex.  On  the  other  hand,  if  the  algorithm  appears  to 
be  struggling,  the  user  is  also  given  the  option  of  restarting 
around  the  best  point  discovered.  That  is,  the  best  point  is  now 
denoted  and  initialization  option  b)  (defined  on  p.  16)  is 

implemented. 
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IV 


PARAMETER  AND  SUBROUTINE  SPECIFICATIONS 


The  (FORTRAN)  program  consists  of  a  main  program,  seven 
general  subroutines,  and  two  user  supplied  subroutines.  In  greater 
detail  these  are  described  next. 

SUBROUTINE  BOX  is  called  by  the  main  program  to  set  up 
the  initial  complex,  call  the  basic  iterative  subroutine  GOBOX,  and 
print  the  final  answer. 

SUBROUTINE  GOBOX  is  the  basic  iterative  subroutine;  it 
accepts  an  initial  complex  and  iteratively  constructs  new  ones 
until  termination.  As  necessary  it  calls  the  special  purpose  sub¬ 
routines  MOVE,  FEASBL,  CENTER,  and  the  function  evaluator  FUNC. 

SUBROUTINE  MOVE  repeatedly  defines  a  new  point  halfway 
between  two  given  points  until  certain  termination  criteria  are 
met. 

SUBROUTINE  FEASBL  checks  that  a  given  vector  belongs  to 
sets  ft  of  the  form  (2) .  If  not  it  attempts  to  reposition  this 
vector  so  as  to  meet  all  constraints. 

SUBROUTINE  CENTER  computes  the  centroid  of  a  complex 
less  one  point. 

SUBROUTINE  GBO  defines  an  all  random  initial  complex 
according  to  initialization  option  c) . 

SUBROUTINE  GDI  defines  a  random  complex  to  include  one 
given  feasible  point  according  to  initialization  option  b) . 

SUBROUTINE  FUNC  computes  the  value  of  a  given  function 


at  a  given  point  in  its  domain. 

SUBROUTINE  IMPLCT  defines  the  various  implicit  constraint 
functions  appearing  in  (2). 

Program  parameters  and  control  variables  are  user  supplied 
as  needed  through  a  NAMELIST  specification  statement.  It  assigns 
the  name  "GO"  to  refer  to  the  following  list  of  variables  and 
array  names... 

Problem  description: 

NX:  number  of  independent  variables 

NIC:  number  of  implicit  constraint  functions  g..  in  (2) 

MAXMIM:  set  «  +1(-1)  to  maximize  (minimize) 

XMIN:  array  of  lower  bounds  for  independent  variables 

XMAX:  array  of  upper  bounds  for  independent  variables 

XIMIN:  array  of  lower  bounds  on  constraint  functions 

XIMAX:  array  of  upper  bounds  on  constraint  functions 

Basic  iteration  controls: 

KPTSET:  used  to  prescribe  number  of  points  in  complex 

N1CUTS:  upper  bound  on  number  of  halfway  moves  towards  centroid 

N2CUTS:  upper  bound  on  number  of  halfway  moves  away  from  centroid 

STEP:  step  size  of  reflection  in  basic  iteration  (see  (3)) 

NLOOPS:  upper  bound  on  number  of  basic  iterations 

convergence  parameter  =  number  of  consecutive  iterations 
without  sufficient  progress. 


NTOL: 


ABSTOL: 


convergence  parameter  (see  (4)  ) 
RELTOL:  convergence  parameter  (see  (5)) 


Initialization  and  run  parameters: 


IRUN: 
IPRINT 
MAX  RAN 
NRANPT 

SEED: 

XO: 

XXO : 


program  operation  control;  usage  defined  in  section  V 

output  control;  usage  defined  in  section  III 

bound  on  number  of  attempts  to  establish  initial  complex 

total  number  of  random  draws  from  feasible  set  to  define 
initial  all-random  complex 

seed  used  in  system  uniform  random  number  generator 
initial  feasible  point 
initial  feasible  complex 


Of  these  variables  it  is  imperative  that  the  user 
assign  values  to  the  first  five  and,  if  NIC  >  0,  to  the  next  two 
also.  Having  done  so  the  program  will  run  using  the  default  set¬ 
tings  of  the  remaining  parameters;  in  particular  an  all-random 
initial  complex  will  be  chosen. 

In  practice,  the  user  will  want  to  consider  the  appro¬ 
priate  values  for  the  initialization  and  convergence  parameters. 
According  to  the  initialization  option  chosen,  one  of  XO,  XXO  or 
NRANPT  must  be  defined.  The  user  will  also  probably  need  to  assign 
relevant  values  to  the  convergence  tolerance  values  ABSTOL  and 
RELTOL,  and  the  iteration  bound  NLOOPS.  By  default,  ABSTOL-O. ,RELTOL= 
0.000001  and  NLOOPS=500.  Also,  by  default,  STEP=1.5, 
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KPTS=IFIX  (1.5*NX+0. 5) ,  N1CUTS=8,  N2CUTS=16,  and  IPRINT=1.  The 
use  of  the  program  control  IRUN  is  described  in  the  next  section, 
along  with  several  examples. 

As  already  noted  the  STEP  parameter  should  satisfy 
1.0  <  STEP  £  2.0.  Ideally  one  would  like  STEP  larger  initially  to 
make  rapid  progress  towards  the  solution  point  and  smaller  as  this 
point  is  approached.  In  general,  larger  values  of  STEP  can  result 
in  infeasible  reflected  points  and  consequent  loss  of  time  in 
retreating  towards  the  centroid,  while  small  values  of  STEP  result 
in  slow  progress  and  increased  iterations.  However,  as  the  length 
of  a  reflection  step  is  roughly  proportional  to  the  diameter  of  the 
associated  complex,  the  step  lengths  automatically  decrease  as  the 
algorithm  progresses  and  the  complexes  shrink.  Hence  algorithm 
performance  is  relatively  indifferent  to  moderate  changes  in  STEP 
and  STEP  =  1.5  has  proved  to  be  a  satisfactory  default  value. 

Finally,  there  is  the  matter  of  "relevant  values"  for  the 
tolerances  ABSTOL  and  RELTOL.  In  general,  only  one  of  these  needs 
be  positive.  If  the  user  has  a  good  idea  of  the  order  of  magnitude 
of  the  optimal  function  value  (as,  for  instance,  in  Example  4  of 
the  next  section) ,  then  ABSTOL  should  be  set  according  to  the  final 
accuracy  desired.  If,  on  the  other  hand,  the  user  has  no  reliable 
guide  to  the  magnitude  of  the  final  function  values,  we  suggest  the 
default  values  of  these  two  tolerance  parameters.  These  cases  are  illustrated 
in  Examples  5,  6  belcw  where,  in  fact,  even  smaller  values  of  RELTOL  are  employed. 


V.  PROGRAM  OPERATION  AND  EXAMPLES 


To  utilize  the  program  for  the  solution  of  a  particular 


optimization  problem  of  the  form  (1) ,  (2) ,  the  user  should  do  the 


following: 

1) 


2) 

3) 


IRUN  =  0 

1 

2 

3 

4 


5: 


6 : 
7: 


define  the  objective  function  f  in  SUBROUTINE  FUNC, 
and  the  constraint  functions  g^  in  SUBROUTINE  IMPLCT. 
In  the  latter  the  notation  XI (J)  is  used  for  g  ^  , 
with  J=l,  2,  . ..,  NIC; 

alter  the  DIMENSION  statement  in  the  main  program, 
if  necessary,  to  accomodate  values  NX>10  or  NIC>9. 
Compile  and  run,  then: 

input  desired  values  of  the  relevant  NAMELIST 
variables,  finally  assigning  a  value  to  the  control 
IRUN  from  among  the  following  options. 

stop 

begin  with  user  specified  complex  XXO 

begin  with  user  specified  feasible  point  XO 

begin  with  all-random  complex  (specify  NRANPT) 

continue  with  further  iterations  from  most  recent 
complex 

restart  as  in  IRUN  =  2  with  XO  set  equal  current 
best  point 

reset  all  NAMELIST  variables  to  default  values 
list  current  value  of  all  NAMELIST  variables 


We  now  illustrate  the  operation  of  the  program  and  the 
results  obtained  on  a  variety  of  examples.  The  actual  outputs 
reported  below  result  from  operation  of  a  double  precision  FORTRAN 
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version  of  our  method  in  CMS  on  an  Amdahl  470/V7  mainframe  computer. 
A  listing  of  this  version  appears  in  the  Appendix. 

2  2  2 

Example  1.  Minimize  100 (x2-x^  )  +  (1-x^)  . 

This  is  a  famous  test  problem  associated  with  the  name  of  Rosen- 
brock.  It  is  an  unconstrained  minimization;  the  objective  function 
surface  is  notable  for  a  narrow  deep  curving  valley.  The  true 
solution  is  clearly  0  at  (1,1).  We  implement  the  program  by 
setting  NX=2,  NIC=0,  XMIN=2*-2„  XMAX=2*2.  (thereby  constraining 

(Xj^Xj)  to  the  cube  defined  by  | x ^  j  <_2 )  ,  and  MAXMIN  =  -1.  We 
report  in  Fig.  4  below  the  results  of  three  sets  of  10  runs  with 
the  indicated  values  of  ABSTOL,  RELTOL,  and  STEP.  (In  all  the 
examples  the  default  value  STEP  =  1.5  is  used  unless  otherwise 
noted.)  In  all  cases  the  conventional  starting  value  XO  =  (-1.2, 
1.0)  was  used  (.IRUN  =  2)  .  For  this  example,  and  in  general  when  we  can 
reasonably  guess  at  the  magnitude  of  the  optimal  function  value,  we  should  set 
RELTOL  =  0.0  and  choose  ABSTOL  according  to  the  desired  accuracy.  We  did  so 
in  the  second  set  of  runs  here,  desiring  that  F  be  reduced  belcw  10  i0.  The 
third  set  is  the  same  except  for  STEP  =  1.3  (the  original  Bax  recarmenda t ion ) ; 
in  this  case  we  note  a  slightly  poorer  and  more  erratic  result. 


Function  Fined. 

Evaluations  Function  Value 


ABSTOL  =  0.0 

RELTOL  =  0.000001 

916 

(323) 

WEEEEWM 

ABSTOL  =  5.0D-11 

360 

RELTOL  =  0.0 

(53) 

ABSTOL  =  5.0D-11 

317 

6.773  D-09 

RELTOL  =  0.0 

STEP  =  1.3 

(63) 

. 

(18.333)D-09 

Fig.  4.  Means  and  (standard  errors)  for  the  Rosenbrock  function  minimization. 
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We  note  that  since  feasible  points  are  readily  available  in 
this  example,  the  pure  random  start,  being  wasteful  of  function 
evaluations,  is  avoided. 

2  2 

Example  2.  Maximixe  x^  +4x^x2+7x2  over  the  set  ft={(x^,x2): 
0<x^<_l,  -1<x^+2x2<_1  ,  -1<_3x^-4x2_<1  } .  This  is  a  concave  quadratic 
program  over  the  pentagonal  region  ft  with  exact  solution  1.48  at 
the  vertex  (0.2,  0.4).  We  set  NX=NIC=2,  XMIN=2*0.1,  XMAX=2*1.1, 
XIMIN=2*-1 . 1 ,  XIMAX=2*1 . 1 ,  ABSTOL=RELTOL=. 000001 ,  and  made  ten 
runs  starting  from  XO=(0.3,  0.2).  The  results  were 


function 

function 

solution 

evaluations 

value 

point 

mean 

54 

1.478883 

(.200004,  .399998) 

(standard  error) 

(10) 

(0.000012) 

((.000010,  .000005)). 

2  4  6  2  4 

Example  3.  Minimize  4x^  -2.1x^  +x^  /3+x^x2~4x2  +4x2  over 
the  rectangle  ft  =  {  (x^,  x2)  :  |XjJ<2.5,  [  x2  |  <_1 . 5  } .  This  is  another 

famous  test  problem,  the  objective  function  being  known  as  the  "six 
hump  camel  back  function"  [17] .  It  is  a  challenge  for  local 
gradient-utilizing  methods  because  of  the  presence  of  an  assortment 
of  stationary  points:  6  local  minima,  2  local  maxima,  and  7  saddle 
points.  Published  numerical  tests  of  other  global  optimization 
procedures  [17]  give  a  global  minimum  value  of  -1.0316  at  the 
points  (0.08984,  -0.71266)  and  (0.71266,  -0.08984),  at  a  cost  of 
72-800  function  evaluations.  We  set  ABSTOL=RELTOL=0 . 0001  and  made 
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5  runs  each  with  KPTS=3,  KPTS=4,  starting  in  all  cases  at  X0=(0,0) 
(IRUN=2) .  The  results  were 


function  evaluations  function  value 


KPTS=3 

132 

-1.031626 

(21) 

(0.000001) 

KPTS=4 

130 

-1.031617 

(14) 

(0.000012) 

We  note  about  the  same  amount  of  work  involved  with  the  larger 
complex,  but  also  a  slightly  greater  variability  in  the  final 
answer . 

Examole  4 .  Solve  the  nonlinear  system  of  equations 

i  3  3 

x2  =  x2 

2 

6xl  “  x2  +  *2  =  °* 

By  sketching  the  corresponding  pair  of  curves  in  the  (x1#x2)  - 

plane,  one  can  easily  verify  a  pair  of  exact  solutions  at  (0,0)  and 

(2,4)  and  one  other  solution  near  (1.5,  -  3.0).  To  determine  this 

third  solution  exactly  we  try  to  minimize  the  sum  of  the  squared 

3  3  2  2  3 

residuals;  that  is,  minimize  (2x^  x2~x2  )  +  (6x^-x2  +x2)  .  The 

least  value  is  of  course  0,  and  the  solution  point  will  be  our 
unknown  third  solution.  We  set  ft={(x^,  x2):  -  4<x2£-2), 
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X0= (1.3,  -  3.0),  and  ABSTOL  =  1.0D-13,  RELTOL  =  0.0,  (attempting  to 
obtain  six  or  seven  significant  figures) .  The  results  of  five 


runs  were 


function  evaluations 

solution  point 

mean 

291 

(1.464352,  -2.506013) 

(standard  error) 

(24) 

(0.000000,  0.000000) 

In  all  cases  the  objective  function  was  reduced  to  less  than 
5.0D-12,  with  an  average  final  value  =  8.85D-13. 

2  2  2 

Example  5.  Minimize  x.^  +x2  +x3  over  the  region 
n={(x.,  x2,  x^)  :  |xi|_<10f  3_<x1x2x3,  3£x1+x2~x3}.  This  problem  is 
unusual  because  fi  is  disconnected:  it  consists  of  three  components 
each  of  equal  distance  from  the  origin.  Hence  there  are  three 
distinct  solution  points,  each  of  which  is  a  boundary  point.  These 
points  can  be  obtained  via  the  Kuhn-Tucker  conditions  and  analysis 
of  the  possible  sign  patterns  of  the  coordinates,  viz.,  (+,+,+), 
(+,-,-),  and  (-,+,-).  Thus  the  first  point  must  obey  the  conditions 


x 


X1  = 
X1X2X3 


1+X2_X3 


x 


2 


3 

3. 


2 

This  leads  to  the  cubic  equations  x  (2x-3)=3  for  x^;  its  solution 
=  1.910820  =  x24a.  Therefore,  x3=2x^-3  =  0. 8216404b.  By  symmetry, 
the  other  solution  points  are  (a,  -b,  -a)  and  (-b,  a,  -a) ,  and  the 


common  minimum  function  value  =  7.977559 


This  example  is  meant  to  typify  the  situation  wherein  we 
are  ignorant  of  the  geometry  of  the  feasible  set  52.  In  such  cases 
we  must  resort  to  the  pure  random  start  option  (IRUN=3) .  Setting 
NRANPT  =  500,  ABSTOL  =  0.,  and  RELTOL  =  1.0D-08,  the  program  easily 
reveals  the  three  distinct  solution  points.  The  results  of  ten 
runs  were 


function  evaluations 

function  value 

mean 

1043 

7.977583 

(standard  error) 

(253) 

(0.000046) 

Eight  of  the  ten  runs  actually  resulted  in  function  values  £7.977560. 

2  2  2  2  2 

Example  6.  Minimize  100  (x^x^  )  +  (1-x^)  +  90(x4~x3  ) 

+  (l-x3)2  +  10.1[(x2-l)2  +  (x4-l)2]  +  19.8 (x2-l) (x4-l) ,  over  the 
four-dimensional  cube  52={  (x1#x2 ,x3,x4)  :  |x^  |£10}.  This  is  another 
test  problem  associated  with  the  name  of  Wood  [6,  p.  403];  it  is 
designed  to  have  a  non-optimal  stationary  point  with  a  corresponding 
function  value  of  approximately  8.0  that  can  cause  premature 
convergence.  The  optimal  solution  is  clearly  0  at  (1,1, 1,1).  A 
starting  value  XO= (-3,-1, -3,-1)  is  suggested.  Our  program  handled 
this  problem  easily  with  no  premature  stalls.  The  results  of  ten 
runs  with  ABSTOL=0.,  RELT0L=1 . 0D-10  (the  value  of  the  objective 
function  was  increased  by  1)  uniformly  yielded  function  values 


<  1.0000000005.  with  an  average  of  1145  (198)  function  values 


Example  7.  Our  final  example  is  a  nonlinear  regression 
taken  from  [24] .  The  problem  is  to  model  the  resistance  (R)  of 
a  thermistor  as  a  function  of  temperature  (T)  via  the  model 


R  =  a  exo 


(tvt)- 


Specifically  we  want  to  assign  values  to  the  parameters  a,b,c 
so  as  to  achieve  the  best  fit  to  the  data 


T 

R 

50 

34,780 

55 

28,610 

60 

23,650 

65 

19,630 

70 

16,370 

75 

13,720 

80 

11,540 

85 

9,744 

90 

8,261 

95 

7,030 

100 

6,005 

105 

5,147 

110 

4,427 

115 

3,820 

120 

3,307 

125 

2,872 
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We  do  so  by  minimizing  the  mean  squared  residual 


The  source  cited  reports  that  this  problem  caused  considerable 
difficulity  to  the  algorithms  being  tested.  In  general  we  note 
that  the  special  structure  of  nonlinear  least  squares  problems  can 
be  utilized  in  the  design  of  specialized  algorithms  such  as  the 
methods  of  Gauss-Newton  or  Levenberg-Marquardt  [25].  For  the 
latest  in  such  algorithms,  see  [26] .  These  may  certainly  be 
expected  to  be  far  more  efficient  than  our  present  all-purpose 
method.  Nevertheless,  it  is  interesting  to  test  our  method  on 
such  a  problem,  recalling  particularly  the  discussion  of  the 
computer  time  vs.  analyst  time  trade-off. 

The  suggested  initial  values  are  XO=(0.2,  4000,  250.), 
the  corresponding  function  value  being  41,153.  We  set  STEP=2. 
and  ABSTOL  =  .0001.  =  RELTOL  .  We  also  introduce  the  artifical 
constraint  g^(a,b,c)  ~  b/(50+c)  £  70  in  order  to  avoid  overflow 
problems  in  the  exponential.  The  outcomes  of  ten  runs  were 
classified  as  success  or  failure  according  as  convergence  was 
or  was  not  achieved  with  at  most  five  restarts  (IRUN=5) .  The 
results  were 


33 
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number 

function 

evaluations 

function 

value 

CPU  time 
(sec. ) 

6480 

9.3788 

B9 

success  6 

■■■ 

(1250) 

( .0010) 

mm 

3031 

267.74 

1.24 

failure  4 

(1340) 

(90) 

(.70) 

The  mean  solution  point  for  the  six  successes  was 

(.005612,  6181.00,  345.210) 

with  standard  errors  (.000006,  0.85,  0.030);  this  may  be  compared 
with  the  optimal  solution  reported  in  [24]: 


(.005609,  6181.,  345.2) 


VI 


CONCLUSIONS 


i 
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% 
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The  numerical  practice  of  nonlinear  multivariate 
optimization  is  an  important,  difficult,  and  much  studied  subject. 
Numerous  algorithms  have  been  proposed,  some  exotic,  some  effective, 
most  requiring  special  hypotheses  on  the  objective  and  constraint 
functions.  We  have  presented  here  an  extremely  simple  and  robust 
procedure  for  handling  the  most  general  nonlinear  inequality  con¬ 
strained  optimization,  and  demonstrated  its  operation  and  effective¬ 
ness  on  a  variety  of  test  problems  of  low  dimension.  The  method 
utilizes  the  basic  Box  complex  iteration  with  modifications  to 
avoid  traps.  The  computer  implementation  provides  a  variety  of 
initialization  and  continuation  features  and  is  designed  for 
interactive  use.  While  problems  of  a  special  nature  can  be  solved 
more  efficiently  with  specially  formulated  algorithms,  in  a  certain 
practical  sense  our  method  is  reasonably  competitive  with  these, 
and  moreover  constitutes  a  safe  reliable  procedure  for  use  on 
problems  with  little  or  no  exploitable  structure.  Examples  of 
the  latter,  on  which  our  method  has  been  successfully  applied, 
include  minimizing  a  discontinuous  cost  function  subject  to  a 
communications  network  reliability  constraint,  and  minimizing  a 
two-stage  survival  probability  function  to  achieve  an  optimal 
attack  against  a  given  layered  ballistic  missile  defense  system. 
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This  appendix  contains  the  complete  FORTRAN  listing  of  the 
optimization  program  described  above.  The  data  for  the  function 
and  implicit  constraint  portions,  SUBROUTINES  FUNC  and  IMPLCT 
respectively,  correspond  to  Example  7  in  Section  V. 


IMPLICIT  REAL*8  <A-Hf  0-2) 

DIMENSION  XX<10f15> fXSTAR(IO) fXMIN< 10) fXMAX< 10) »X0(10) 

DIMENSION  XIMIN<9)fXIMAX<9) fXI<9) fFF<15) fXX0<10f15>  fXCIO) 

COMMON  /BOXCOM/  STEP  t RELTOL f ABSTOL  f NTOL f  N1CUTS  r N2CUTS  ? NLOOPS 

*  fIPRINTfMAXMINfSEEDfNRANPTfMAXRAN 
REAL*8  LXMIN  f  LXMAX f LX IMIN f LX I MAX f LXO f LXXO f LBLANK 
NAMELIST  /GO/  NXfNICfMAXMINfXMINfXMAXfXIMINfXIMAXf 

*  NPTSET f  N1 CUTS  f  N2CUTS f  STEP  f  NLOOPS  f  NTOL  f  ABSTOL  f  RELTOL  f 

*  IRUN  f I PR INT  f  NRANPT f  MAXRAN  f  SEED  f  XO  f  XXO 

DATA  LXMINfLXMAX  /8H  XMIN  f8H  XMAX  / 

DATA  LXIMINfLXIMAX  /8H  XIMIN  f8H  XIMAX  / 

DATA  LXO f LXXO f LBLANK  /8H  XO  f8H  XXO  f8H  / 

C 

C  PROGRAM  NEWBOX  USES  SUBROUTINE  BOX  TO  SOLVE  THE  PROBLEM 
C  MAXIMIZE  MAXMIN#FUNC  <  X ) 

C  SUBJECT  TO  XMIN< I )  .LE.  X(I>  .LE.  XMAX(I)  fI=1f...fNX 

C  XIMIN(J)  .LE.  XI(J)  .LE.  XIMAX(J)  fJ=1f...fNIC 

C  WHERE  XI (J)  IS  THE  VALUE  OF  THE  JTH  IMPLICIT  CONSTRAINT  FUNCTION 
C  EVALUATED  AT  THE  POINT  Xf  AS  COMPUTED  IN  SUBROUTINE  IMPLCT. 

C 

C  THE  USER  MUST  SUPPLY  SUBROUTINE  FUNC(XfNXfF)  WHICH  CALCULATES  Ff 
C  THE  FUNCTIONAL  VALUE  OF  A  POINT  Xf  AND  A  SUBROUTINE 
C  IMPLCT(XfNXfXIfNIC)  WHICH  CALCULATES  XI f  THE  NIC  IMPLICIT 
C  CONSTRAINT  VALUES  CORRESPONDING  TO  THE  POINT  X. 

C 

C 

C  SET  PARAMETERS  TO  THEIR  DEFAULT  VALUES. 

C 

100  KPTSET*0 
N1CUTS=8 
N2CUTS*16 
STEP* 1.5 
NL00PS=500 
NT0L.=5 
ABSTQL=0 . 

RELTOL*. 000001 
I RUN* 3 
IPRINT*1 
NRANPT*0 
MAXRAN=100 
SEED* 184287807 
C 
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C  READ  NAMELIST  &G0  DATA  DEFINING  PARAMETERS  FOR  NEXT  RUN 

C 

200  WRITE <4  f  9  > 

9  FORMAT ( ' 1ENTER  NAMELIST  SGO  INPUT  <IRUN=7  LISTS  ALL)') 
READ ( 5 1 GO  > 

IF  <  I  RUN  ♦  l.E  .  0 )  GO  TO  1000 
IF ( I RUN- 6 )  300 » 100 f 400 
C 

C  CALL  BOX  TO  OPTIMIZE  FOR  THE  CURRENT  PARAMETERS 
C 

300  KPTS-NX+ (NX+1 )/2 

IF ( KPTSET . GE . 2 >  KPTS-KPTSET 

CALL  BOX<XXfNXfKPTSfNICfXMINfXMAXfXIMINfXIMAX 
*  r  FF f  X I f  XSTAR  f  FSTAR  fXfXOfXXOf I RUN ) 

GO  TO  200 
C 

C  PRINT  THE  CURRENT  VALUES  OF  THE  PARAMETERS 


C 

400  WRITE ( 4  f  1 )  NX f  NIC  f MAXMIN 

1  FORMAT ( 'OCURRENT  VALUES  OF  PARAMETERS  ARE ?  '  f  /  f 


*  'O' fIOXf 'PROBLEM  DESCRIPTION '*/ f 

*  '  NXfNICfMAXMIN  =  ' f3I5> 

WR I TE < 6 » 2 >  LXMIN» (XMIN(I) f I  =  lfNX) 


WRITE  <  6  » 2 )  LXMAX  t ( XMAX  <I)rI  =  l*NX) 

WRITE (6. 2)  LXIMIN»<XIMIN<I)f I=1»NIC) 

WRITE <4 »2)  LX I MAX » ( XI MAX (I) rI  =  lrNIC) 

2  FORMAT (A8f /=  '  » 1P10E12 . 4 > / f < 10X > 1P10E12 . 4) ) 

WR I TE ( 6 , 3 )  KPTSET  f N1 CUTS  f  N2CUTS f STEP f NLOOPS  f  NTOL  f ABSTOL  f  RELTOL 

3  FORMAT < 'O' fIOXf 'BOX  PARAMETERS ' f / f 

*  '  KPTSET f N1CUTS f N2CUTS f STEP  =  ' f3I5fF12.6f/f 

*  '  NLOOPSfNTOLfABSTOLfRELTOL  =  ' f2I5f1P2E12.3) 

WR I TE < 6 f  4 )  I RUN  f I FRI NT  f NRANPT  f  MAXR AN  f  SEED 

4  FORMAT ( ' 0 ' f 10X f " RUN  INITIALIZATION  PARAMETERS ' f / f 

*  'IRUNfIPRINTfNRANPTfMAXRANfSEED  *  ' f2I4f2I7f3XfZ16) 
WRITE (6f 2)  LXOf (X0<I) fI=1fNX) 

WRITE <6f 2)  LXXOf (XXO(IfI) fI=1fNX) 

DO  440  J=2 f KPTS 

440  WRITE (4 f 2 )  LBLANKf (XX0(IfJ)fI  =  1fNX) 

GO  TO  200 
C 

1000  STOP 


END 
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SUBROUTINE  BOX(XX,NX,KPTS»NIC»XMIN,XMAX»XIMIN»XIMAX,FF» XI 

*  ,XSTAR»FSTAR,X»XO,XXO» IRON) 

IMPLICIT  REAL*8  (A-H,  0-Z) 

DIMENSION  XX (NX, KPTS) »XMIN(NX>  »XMAX(NX) .XIMIN(NIC) »XIMAX(NIC) 
DIMENSION  FF(KPTS) , XI (NIC) ,XSTAR(NX) 

DIMENSION  X(NX) *XO(NX),XXO(NX,KPTS) 

COMMON  /BOXCOM/  STEP » RELTOL  , ABSTOl.  ,  NTOL  » N1 CUTS  »  N2CUTS , NLOOPS 

*  » I PR I NT  ,HAXM IN, SEED, NRANPT 
REALMS  DDATE , DTTME 

BOX  SETS  UP  THE  INITIAL  COMPLEX  THEN  CALLS  OOBOX  TO  FIND  THE 
OPTIMAL  NX-VECTOR  XSTAR  AND  ITS  FUNCTIONAL  VALUE  FSTAR. 

GO  TO  <100*200»300»400»500) , I RUN 

IRUN-1  ....  SET  INITIAL  XX  TO  USER-DEFINFD  XXO 

100  WRITE ( A » 1  )  KPTS  » STEP 

1  FORMAT < 'CENTERING  BOX  WITH  DEFAULT  XX',/ 

*  KPTS »  STEP  *  ' » I5,F10. A) 

DO  150  ,1=1  » KPTS 

DO  140  1=1 *NX 
140  XX(I,J)=XXO(I, J) 

CALI.  FUNC(XX(1,J)»NX,FF(J)) 

150  CONTINUE 
NFIJNC  =  KPTS 
GO  TO  900 


IRUN=?  ....  SET  INITIAL  XX  RANDOMLY  AROUND  (FEASIBLE)  XO 

200  WRITE ( 6 » 2 )  SEED » KPTS > STEP , XO 

2  FORMAT< CENTERING  BOX  WITH  RANDOM  XX,  SEED  «  ',D13.5,/ 

*  KPTS, STEP  »  ' ,15, FI  0.6,/ 

*  ,'  XO  =',1P?E12.4,/,C9X,1P9E12.4>> 

CALL  GB1 <XX»NX,KPTS,XO,FF,XMIN»XMAX,XIMIN»XIMAX,XI 

*  , NIC, IOK  ) 

NFUNC*KPTS 

IF(IOK.GT.O)  00  TO  900 
GO  TO  1000 

IRUN=3  ....  SET  INITIAL  XX  RANDOMLY  IN  FEASIBLE  RFOION 

300  WRITE  <  A , 3 )  REED, KPTS, NRANPT, STEP 

3  FORMAT< 'CENTERING  BOX  WITH  ALL  RANDOM  XX,  SEED  »  ',D13.5,/ 

*  »'  KPTS, NRANPT, STEP  ■  '»2I5,F10.A> 

CALL  GBO(XX,NX,KPTS,FF,XMIN,XMAX, XIMIN, XIMAX 

*  , XI, NIC, X, IOK) 

NFUNC=MAXO< KPTS, NRANPT) 

IF(IOK.OT.O)  00  TO  900 
00  TO  1000 
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IRUN  =  4  ....  CONTINUF  WITH  XX  AND  FF  AS  I.FFT  BY  PREVIOUS  RUN 

400  WRITE ( A » 4  )  KPTS t NLOOPS > STFP 

4  FORMAT ( 'OCONTINUING  BOX  WITH  SAME  XX' >/ 

*  KPTS » NI.OOPS  *  STEP  =  'r2I5»F10.A> 

NFUN0=0 

GO  TO  900 

IRUN=5  ....  SET  INITIAL  XX  RANDOMLY  AROUND  LAST  XSTAR 

500  WRITF (6,6)  SEED . KPTS , STEP » XSTAR 

6  FORMAT! 'ORESTARTING  BOX  AROUND  XSTAR>  SEED  ='»Dls.5»/ 

*  >'  KPTS » STEP  ='»I5rF10.A»/ 

*  XSTAR  =' »lP9E12.4»/» ( 17X»1P9E12.4) ) 

CALL  GB1 <XXt NX »KPTS» XSTAR? FF ? XMIN ? XMAX ? X JMIN ? XIMAX ? XI 

*  ?  NIC » I OK ) 

NFUNC=KPTS 

IF(IOK.GT.O)  GO  TO  900 
GO  TO  1000 

NOW  CALL  BOX  WITH  XX  AND  FF  TO  FIND  OPTIMUM  XSTAR  AND  FSTAR 

900  CALL  TIMES(PDATE»DTIMFt  IVCPUJ  ?  ITCPIJ) 

CALL  GOBOX(XX»NXfKPTSfNIC?XMINrXMAX?XIMIN?XIMAX?FF?XI 

*  , XSTAR » FSTAR ? MLOOPS ? NFUNC > 

FSTAR=MAXMIN*FSTAR 

CALL  TIMES ( DDATE »  DTIME » IVCPU2? ITCPU) 

TIME3 » 00001 IVCPU2-IVCPU1  ) 

IF ( IPRJNT .GT .0)  WRI TE ( A* 5 >  MLOOPS? NFUNC? TIME? FSTAR? XSTAR 

5  FORMAT ( '0 ' ? 15?  '  LOOPS? '?IA?'  FUNCS?  AND'?F7,3?'  SECS  LATER  BOX' 
*»'  FINDS  FSTAR  AND  XSTAR  OF ' ? 1PF1 A . 9 ? / ? ( 1P9E1 4 . A ) ) 

1000  RETURN 
END 
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SUBROUTINE  GOBOX < XX  *  NX » KPTS *  NIC » XMIN* XMAX » XIHIN * XIMAX 

*  *  FF »  X I »XC»FSTAR»  NLOOPS  r NFUNC ) 

IMPLICIT  REAL  *8  <A-H»  O-Z) 

DIMENSION  XX<NX» KPTS) » XMIN(NX) »XMAX(NX) »XIMIN<NIC) .XIMAX(NIC) 
DIMENSION  FF<KPTS)*XI(NIC)» XC<NX ) 

COMMON  /BOXCOM/  STEP  *  REL  TOI. .  ABSTOL  *  NTOL .  N1CUTS*  N2CUTS*  NLOOPS 

*  t IPRINT  rMAXHTN 
DATA  IO»II  /  0».t  / 

GIVEN  AN  INITIAL  COMPLEX  OF  KPTS  NX-VECTORS  IN  XXt  AND  KPTS 
FUNCTION  EVALUATIONS  IN  FF »  GOBOX  WILL  TRY  TO  FIND  THE  POINT 
WHICH  MAXIMIZES  FUNC  ON  THE  FEASIBLE  REGION.  THIS  POINT  IS 
RETURNED  IN  ARRAY  XC*  AND  ITS  VALUE  IN  FSTAR . 

FIND  HI  AND  1.0  POINTS  IN  THE  COMPLEX 


JLO*  1 
JHT  =  1 

DO  40  J*2 »  KPTS 

IF<FF<  J)  .LT.FF1  Jl.O>  >  JLO-J 
IF(FF(J) . GT .FF  <  JHI >  >  JHI=J 
40  CONTINUE 
I  TOI. -0 

IF  < IPRINT-2)  X00,50t60 
50  F=MAXMIN*FF< JHI) 

WRITE  <6»2)  IO»Il»JHI»F»(XX(ItJHI>»I*i»NX) 

GO  TO  100 
60  DO  70  J®1»KPTS 

F=MAXMIN*FF< J) 

IF(J.NE.JHI)  WRITE(A*1)  IO»IltJ»F»(XX(I»J)»I®ltNX) 
IF(J.EG.JHI)  WRITE  (6*2)  10* II  * JfFf (XX( I » J) f I®1 *NX > 

1  FORMAT ( '  LPfNFfJ»F»X='»3I3»lPE12.4»lXf9EJ1.3»/ 

*  f(37Xf1P9EU,3>> 

2  FORMAT <  '  **LP»NF» J»FrX®'  *  313 r 1 PEI 2 . 4 f 1 X » 9E1 1 .3*/ 

*  f <37Xf 1P9E11 .3) ) 

70  CONTINUE 

MAKE  UP  TO  NLOOPS  ITERATIONS*  FACH  TIME  IMPROVING  THE  LOWEST  VALUE 

100  DO  500  I LOOP® 1 » NLOOPS 

LOCATE  NEXT  LOWEST  POINT 

NORMLC®! 

IFUNC-NFUNC 
NEXTl.O®  JHI 
DO  140  .J®  1  f  KPTS 

IF  <  J  ♦  NE  ♦  JLO  .AND.  FF  <  J ).  I.  T .  FF  < NEXTl.O )  >  NEXTLO-J 
140  CONTINUE 

C 
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FIND  CENTROID  AND  CHECK  ITS  FEASIBILITY < 

CALL  CENTER(XX,NX»KPTS» JLO»XC) 

CALL  FEAS8L ( XC » NX, XC*  XMIN,  XHAX » XININ, XIMAX , XI , NIC 

*  » NMOVES » 0 ) 

IF<NMOVES.GT.O)  GO  TO  200 

FOR  FEASIBLE  CENTROID,  LOOK  FOR  NEW  POINT  BY  MOVING  FROM  LO  POINT 
THRU  CENTER  BY  FACTOR  OF  STEP,  THEN  MOVE  TOWARDS  CENTER  IF  NECESSARY 
TO  GET  A  POINT  BOTH  FEASIBLE  AND  BETTER  THAN  NEXTI.O  POINT. 

DO  150  1=1 » NX 

150  XX(I, JLO)=XX(I» JL0)+<1 .+STEP)*<XC< I )-XX(I. JLO) ) 

CALI.  MOVE (  XX  ( 1 » JLO) , NX » XC , NMOVES » NLCUTS , FF <  JLO)  »FF(NEXTLO) 

*  »NFUNC*XMIN*XMAX»XIMIN»XIMAX»XI»NIC) 

IF ( NMOVES .LE. NIC UTS)  00  TO  400 

IF  ABOVE  DIDN'T  WORK,  OR  IF  CENTER  NOT  FEASIBLE,  LOOK  FOR  NEW 
POINT  BY  MOVING  FROM  CENTER  TOWARDS  HI  POINT. 

200  DO  250  1=1, NX 

250  XX  < I , JLO ) =XC ( I ) 

NORMLC=- 1 

CALI.  MOVE  (  XX  <1 ,  JLO ) ,  NX ,  XX  (  1 ,  JHI  )  ,  NMOVES  ,  N2CUTS 

*  ,FF<Jl .0)  pFF (NEXTI.O)  'NFUNC 

*  ,XMIN,XMAX,XIMIN, XIMAX, XI, NIC) 

IF ( NMOVES .LE.N2CUTS)  00  TO  400 

IF  THIS  DIDN'T  WORK,  TEST  ON  BOTH  SIDES  OF  HI  POINT  ALONG  LINE 
THRU  CENTER  FOR  FEASIBILITY  AND  IMPROVEMENT.  IF  BOTH  ARE  INFEAS¬ 
IBLE  OUIT  WITH  INFEASIBLE  DIRECTION  AT  B00»  IF  ONE  IS  BOTH 
FEASIBLE  AND  BETTER,  CONTINUE  WITH  ITI  OTHERWISE  QUIT  WITH  NO 
IMPROVEMENT  AT  700 

CALI.  FEASBl.  <  XX  <  1 ,  Jl.Q ) ,  NX  ,  XX  <  1  ,  JL  O ) ,  XMIN ,  XMAX  »XIMIN»  XIMAX 

*  , XI , NIC, NMOVES, 0) 

IF  <  NMOVES  .  l.E  .  0 )  GO  TO  700 

DO  ?60  1=1, NX 

260  XX ( I >  JLO ) =XX  < I » JHI )  +  ( XX< I , JHI ) -XX< I » JLO) ) 

CALL  FEASBL <XX( 1 , JLO) ,NX , XX ( 1 , JLO) , XMIN, XMAX , XININ, XIMAX 

*  , XI, NIC, NMOVES, 0) 

IF(NMOVES.OT.O)  00  TO  800 

CALL  FUNC ( XX < 1 , JLO ) ,NX » FF ( JLO) ) 

NFUNC=NFUNC+1 

IF(FF( JLO) .OT.FF(NEXTLO) )  GO  TO  400 
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C  NOW  THAT  WE  HAVE  IMPROVED  THE  LOU  X»  TEST  FOR  COMPLETION  AND 
C  LOOP  AROUND  AGAIN 
C 

400  MFUNC*NORMI.C#  <  NFUNC-IFUNC  ) 

F=WAXMIW*FF  ( JLO) 

IF  ( FF  <  JLO )  . I.E  ♦  FF<  JHI )  )  GO  TO  440 

IF  < IPRINT . GE • 2 )  URITE<6»2)  I  LOOP  »  MFUNC » Jl.O » F 

*  *  (XX(I»  Jl.O)  »I  =  lfNX) 

JHI= JLO 

GO  TO  450 

440  IF(IPRINT.GE.3)  URITE(6»1)  IL OOP » MFUNC » JLO » F 

*  »<XX<I»JL0)»I=1»NX> 

450  Jl.O=N£XTLO 

IT0L»IT0L+1 

ATOI..*DABS <  FF  ( JHI  ) -FF  (  JL  0  > ) 

IF < ATOl. ♦  LE . ABSTOL )  GO  TO  480 
IF  <  ATOl.  •  GT « RELTOL  *DABS <FF< JHJ  > ) >  IT0L*0 
480  IF< ITOL.GE.NTOL)  GO  TO  600 

500  CONTINUE 


MLOOPS=NLOOPS+1 
GO  TO  900 
C 

600  MLOOPS=Il.OOP 
GO  TO  900 
C 

700  Ml..  OOPS*- 1  LOOP 
GO  TO  900 
C 

800  MLOOPS  =  -Nl.OOPS-I 
C 

900  DO  950  1*1 »NX 
950  XC<I)*XX(I» JHI) 
FSTAR*FF< JHI ) 

C 

RETURN 

END 


NO  CONVERGENCE 

CONVERGENCE 

UNABLE  TO  IMPROVE 

INFEASIBLE  DIRECTION 
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SUBROUTINE  FEASBL(X.NX»XC.XMIN» XMAX .XIMIN.XIMAX.XI.NIC 

*  . NMOVES. MAXMOV) 

IMPLICIT  REAL*8  (A-H.O-Z) 

DIMENSION  X(NX)»XC(NX)»XMIN( NX ) .XMAX(NX) 

DIMENSION  XIMIN<NIO) .XIMAX(NIC) » XI (NIC) 

COMMON  /BOXCOM/  STEP »  REl.TOL  .  ABSTOL  .  NTOL » N1  CUTS » N2CUTS » NLOOPS 

*  . IPRINT.MAXMIN 

c 

C  TRIES  TO  MAKE  VECTOR  X  FEASIBLE  BY  (IF  NECESSARY)  MOVING  X 
C  TOWARDS  GIVEN  FEASIBLE  VECTOR  XC  UP  TO  MAXMOV  TIMES. 

C 

C  CHECK  AGAINST  EXPLICIT  CONSTRAINTS 
C 

NM0VES=0 

100  DO  200  1  =  1 » NX 

IF(X(I) .LT.XMIN(I)  .OR.  X ( I ) . GT . XMAX ( I ) )  GO  TO  400 
200  CONTINUE 

CHECK  AGAINST  IMPLICIT  CONSTRAINTS 

IF  (  NIC  .  l.E  .  0  )  RETURN 
CALL  I MPLCT (X.NX.XI.NTC) 

DO  .TOO  1  =  1  .NIC 

IF(XI(I).LT.XIMIN(I)  .OR.  XI ( I ) . GT . XIMAX ( I ) )  GO  TO  400 
300  CONTINUE 

RETURN 

INCREMENT  NMOVES  AND  IF  .LE.  MAXMOV.  MOVE  X  TOWARDS  CENTROID. 

400  NMOVES=NMOVES+ 1 

IF(NMOVES.GT. MAXMOV)  RETURN 
DO  500  1=1. NX 
500  X(I)=0.5*(X(I)+XC(I)) 

GO  TO  100 
C 

END 
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SUBROUTINE  MOVE ( X . NX » XT  » NMOVES >  MAXMOV .  F » EMIN » NFUNC 

*  » XMIN » XMAX » X IMIN » XIMAX » X I »NIC) 

IMPLICIT  REAL  *8  <A-H»0-Z> 

DIMENSION  X ( NX ) »XT<NX) » XMIN(NX) »XNAX< NX >»XIMIN(NIC) 
DIMENSION  XIMAX(NIC) »XI(NIC) 

TRIES  TO  VERIFY  THAT  VECTOR  X  IS  BOTH  FEASIBLE  AND  HAS  VALUE 
.GT.  FMIN.  IF  NECESSARY  BY  MOVING  X  TOWARDS  THE  VECTOR  XT  UP 
TO  MAXMOV  TIMES. 

NM0VES*0 

100  CALL  FEASBl.<X'NX'XT»XMIN»XHAXtXIMIN'XIMAX'XI'NIC 

*  f IMOVES.MAXMOM-NMOVES) 

NMOVES=NMOVES+IMOVES 
IF(NMOVES.GT. MAXMOV)  RETURN 

CALL  FUNC  <  X » NX » F ) 

NFUNC»NFUNC+1 
IF(F.GT.FMIN)  RETURN 
DO  200  1  =  1 »  NX 
200  X<I)»0.5*(X(I)+XT(I)) 

NM0VES=NM0VFS+1 

IF<NMOVES.LE. MAXMOV)  GO  TO  100 
END 


SUBROUTINE  CENTER  <  XX  *  NX  » KPTS » NOTJ » XC ) 

IMPLICIT  REAL*8  <A-HrO-Z> 

DIMENSION  XX < NX f KPTS) *XC<NX> 

COMPUTES  THE  NX-VECTOR  XC  AS  THE  CENTROID  OF  THE  KPTS  NX-VECTORS 
IN  XX  LESS  THE  ONE  INDEXED  BY  NOTJ. 

DIVKM1*1.0/<KPTS-1) 

DO  200  1=1 .NX 
XSUM=0. 

DO  100  J=1 »KPTS 
100  XSUM=XSUM+XX  < I » J ) 

XC<I)=<XSUM-XX<I»N0TJ)>#DIVKM1 
200  CONTINUE 

RETURN 

END 
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SUBROUTINE  GBO (XXfNXf NPTS fFFfXMINf XMAX f XIMIN f XIMAX 

*  f  XI f  NIC  f X  t IOK ) 

IMPLICIT  REAL*8  (A-HfO-Z) 

DIMENSION  XX ( NX f KPTS) fFF (KPTS) fXMIN(NX) f XMAX (NX) 

DIMENSION  XIMIN (NIC) f XIMAX (NIC) fXI (NIC) fX(NX) 

COMMON  /BOXCOM/  STEPfRELT0LfABST0LfNT0LfN1CUTSfN2CUTSfNL00FS 

*  f IPRINT  f MAXMIN  f SEED , NRANPT f  MAXRAN 

C 

C  CREATES  A  RANDOM  INITIAL  COMPLEX f  XX f  BY  CHOOSING  THE  BEST  KPTS 
POINTS  OUT  OF  A  COLLECTION  OF  MAX(KPTSfNRANPT)  RANDOM  VECTORSf 
EACH  UNIFORMLY  DISTRIBUTED  OVER  THE  FEASIBLE  REGION. 

FILL  XX  WITH  KPTS  FEASIBLE  RANDOM  VECTORS 

NRANX=0 

MXRANX=MAXRAN*MAXO (KPTS f  NRANPT ) 

DO  300  J-l f KPTS 
100  DO  200  1=1 fNX 

200  XX ( I f  J ) =XMIN ( I ) +RAN2 ( SEED ) * ( XMAX ( I ) ~XM IN ( I ) > 

NRANX=NRANX+1 

IF ( NRANX . GT . MXRANX )  GO  TO  1100 

CALL  FEASBL(XX(1fJ) fNXfXX(1f J> fXMINfXMAXfXIMINfXIMAXfXI 

*  fNICfNMOVESfO) 

IF(NMOVES.GT.O)  GO  TO  100 
CALL  FUNC(XX(1f J) fNXfFF( J) ) 

300  CONTINUE 

K1=KPTS+1 

IF (Kl.GT. NRANPT)  GO  TO  1000 
JMINF=1 

DO  500  J=2fKPTS 

500  IF ( FF ( J ) . LT . FF ( JMINF ) )  JMINF=J 


C  CHOSE  (NRANPT-KPTS)  RANDOM  VECTORS*  ALWAYS  KEEPING  IN  XX 
C  THE  BEST  KPTS  VECTORS  LOOKED  AT  SO  FAR. 

C 

DO  900  K=K1 »NRANPT 
600  DO  650  1=1 »NX 

650  X< I )=XMIN( I ) +RAN2 ( SEED ) * < XMAX < I )-XMIN< I ) ) 

NRANX=NRANX+ 1 

IF (NRANX.GT .MXRANX)  GO  TO  1100 

CALL  FEASBL ( X » NX r X r XMIN*  XMAX *  XIMIN*  XIMAX *  XI  *  NIC 

*  *NMQVES  1 0 ) 

IF ( NMOVES . GT . 0 )  GO  TO  600 
CALL  FUNC ( X »NX »F > 

IF(F.LE.FF( JMINF) )  GO  TO  900 
FF ( JMINF ) =F 
DO  700  1=1 »NX 
700  XX(If JMINF)=X(I) 

JMINF=1 

DO  800  J=2*KF‘TS 

800  IF(FF( J) .LT.FF( JMINF) )  JMINF=J 

900  CONTINUE 

C 

1000  I0K=1 
RETURN 
C 

1100  WRITE ( 6 f 3 )  NRANX 

3  FORMAT < '0#****  IN  GBO»  WE  HAVE  EXCEEDED  THE  MAX  NUMBER' 

*  f'  OF  RANDOM  POINTS  <'rI6»'  )') 

I0K=0 

RETURN 

END 
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SUBROUTINE  GB1 ( XX » NX , KPTS » XO, FF » XM IN, XMAX , XIMIN, XIMAX ,  XI 

*  ,NIC,IOK) 

IMPLICIT  REAL *8  (A-H*0-Z> 

DIMENSION  XX (NX, KPTS) ,XO(NX> ,FF (KPTS) »XMIN(NX) , XMAX (NX) 
DIMENSION  XIMIN( NIC ), XIMAX (NIC) , XI (NIC) 

COMMON  /BOXCOM/  STEP,RELT0L,ABST0L,NTQL,N1CUTS,N2CUTS,NL00PS 

*  , IPR INT , MAXMIN, SEED » NRANPT , MAXRAN 

CREATES  A  RANDOM  INITIAL  COMPLEX ,  XX,  AROUND  AND  INCLUDING 
A  GIVEN  FEASIBLE  VECTOR  XO. 

CHECK  XO  FOR  FEASIBILITY  AND  THEN  PUT  IT  IN  XX 

CALL  FEASBL  ( XO » NX , XO , XMIN, XMAX , XIMIN, XIMAX  r  XI f NIC 

*  » NMOVES » 0 ) 

IF ( NMOVES .  GT .0)  GO  TO  1000 
DO  50  1=1 , NX 

50  XX(I,1)=X0(I) 

CALL  FUNC(XO»NX»FF (1 )  ) 

C 

C  FILL  REST  OF  XX  WITH  RANDOM  VECTORS  WHICH  HAVE  BEEN  FORCED 
C  FEASIBLE  BY  MOVING  THEM  TOWARDS  XO  IF  NECESSARY. 

C 

NRANX=0 

MXRANX=MAXRAN* ( KPTS- 1 ) 

DO  400  J=2  , KPTS 
100  DO  200  1=1, NX 

200  XX ( I , J ) =XMIN ( I ) +RAN2( SEED ) *  <  XMAX ( I ) -XMIN ( I > ) 

NRANX=NRANX+1 

IF < NRANX . GT . MXRANX )  GO  TO  1100 

CALL  FEASBL<XX< 1 » J) ,NX,X0, XMIN, XMAX, XIMIN, XIMAX, XI » NIC 

*  » NMOVES, N2CUTS) 

IF (NMOVES. GT .N2CUTS)  GO  TO  100 
CALL  FUNC(XX( 1 » J) »NX,FF ( J) ) 

400  CONTINUE 

C 

IQK  =  1 
RETURN 
C 

1100  WRITE (6,2)  NRANX 

2  FORMAT ( '0**#**  IN  GB1 ,  WE  HAVE  EXCEEDED  THE  MAX  NUMBER' 

*  ,  '  OF  RANDOM  POINTS  (',16,'  )') 

I0K=0 

RETURN 

C 

1000  WRITE (6,1)  XO 

1  FORMAT ( '0*****  NON-FEASBIBLE  INITIAL  POINT  GIVEN  TO  6B1 ' ,/ 

*  ,'  XO  =' ,1P10E12.4) 

I0K=0 

RETURN 

END 
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SUBROUTINE  FUNC(X»N»F> 

IMPLICIT  REAL. $8  (A-HfO-Z) 

DIMENSION  X ( N ) 

COMMON  /BOXCOM/  STEP .  REL  TOI.  t  ABSTOI.  *  NTOL  »  N1CUTS  *  N2CUTS » NLOOPS 

*  » IPRINTfMAXMIN 
DIMENSION  Y<16> 

DATA  Y/34780. » 28610. .23650. » 1 9630. » 16370. » 13720. .11540. .9744. »8261 

*  . > 7030. . 6005. » 5147. .4427. . 3820. . 3307. f 2872,/ 

RESO(A.B.C»Y.T>  *  <Y  -  A*DEXP< B/< T+C > > >**2 

C 

C  COMPUTES  F»  THE  VALUE  OF  THE  FUNCTION  WHEN  EVALUATED 
C  AT  THE  POINT  X  (AN  N-VECTOR) . 

C 

F  =  0. 

DO  100  J*l»16 
T  =  50.+5.*( J-l) 

100  F  =  F  +  RESQ(X<1)»X<2)»X<3>»Y(J).T) 

F  =  DSQRT(F) 

C 

F=MAXMIN*F 

RETURN 

END 


SUBROUTINE  IMPLCT (X.NX.XI.NIC) 

IMPLICIT  REAL*8  (A-H.O-Z) 

DIMENSION  X ( NX )»XI(NIC) 

COMPUTES  XI.  THE  NIC-VECTOR  OF  IMPLICIT  CONSTRAINT 
FUNCTION  VALUES  AT  THE  NX-VECTOR  X. 

XI ( 1 )  =  X ( 2 ) / <50. +X  <  3 ) ) 

C 

RETURN 

END 
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